#########################################################################################
# 
# Alternative Mechanism: Post-traumatic Growth
#
#########################################################################################

#########################################################################################
### Democratic Republic of Congo

alpha(cbind(congo$ptg1, congo$ptg2, congo$ptg3, congo$ptg4, congo$ptg5, congo$ptg6, congo$ptg7, congo$ptg8, congo$ptg9, congo$ptg10))

# Mediator Equation
drc.sv.ptg <- ictreg.joint(Y ~ female + age.z + edu.z + income.z  + hh_size.z +  murder_yes 
                         + terr_2 + terr_3 + terr_4 + terr_5 + terr_6, 
                         J=3, 
                         data=congo[is.na(congo$activity_prev)==F, ], 
                         treat="treatment",
                         outcome="PTG",
                         outcome.reg="linear",
                         constrained=TRUE,
                         maxIter = 20000)
summary(drc.sv.ptg)

# Outcome Equation
drc.ptg.sp <- ictreg.joint(Y ~ female + age.z + edu.z + income.z  + hh_size.z +  murder_yes 
                         + terr_2 + terr_3 + terr_4 + terr_5 + terr_6 + activity_prev + PTG, 
                         J=3, 
                         data=congo, 
                         treat="treatment",
                         outcome="soc_part.2",
                         outcome.reg="linear",
                         constrained=TRUE,
                         maxIter = 30000)
summary(drc.ptg.sp)

# Predicted Values
predict.ptg <- predict.ictreg.joint(drc.sv.ptg, sensitive.value="both", return.draws=T)

X.1 <- as.data.frame(drc.ptg.sp$x)
X.1$PTG <- predict.ptg[[1]]$Z1[,1]

X.0 <- as.data.frame(drc.ptg.sp$x)
X.0$PTG <- predict.ptg[[1]]$Z0[,1]

predict.sp.1 <- predict.ictreg.joint(drc.ptg.sp, newdata=X.1, sensitive.value="both", return.draws=T)
predict.sp.0 <- predict.ictreg.joint(drc.ptg.sp, newdata=X.0, sensitive.value="both", return.draws=T)

# Y(M=1, T=1)
Y.11.drc <- mean(predict.sp.1[[1]]$Z1[,1])

# Y(M=1, T=0)
Y.10.drc <- mean(predict.sp.1[[1]]$Z0[,1])

# Y(M=0, T=1)
Y.01.drc <- mean(predict.sp.0[[1]]$Z1[,1])

# Y(M=0, T=0)
Y.00.drc <- mean(predict.sp.0[[1]]$Z0[,1])

# ACME
Y.11.drc-Y.01.drc
Y.10.drc-Y.00.drc
((Y.11.drc-Y.01.drc)+(Y.10.drc-Y.00.drc))/2

# ANDE
Y.11.drc-Y.10.drc
Y.01.drc-Y.00.drc
((Y.11.drc-Y.10.drc)+(Y.01.drc-Y.00.drc))/2

# With Uncertainties
# Y(M(1), T=1)
Y.11.drc <- apply(predict.sp.1$draws.predict[[2]], 2, mean)
mean(Y.11.drc)
quantile(Y.11.drc, c(.025, .975))

# Y(M(1), T=0)
Y.10.drc <-  apply(predict.sp.1$draws.predict[[1]], 2, mean)
mean(Y.10.drc)
quantile(Y.10.drc, c(.025, .975))

# Y(M(0), T=1)
Y.01.drc <- apply(predict.sp.0$draws.predict[[2]], 2, mean)
mean(Y.01.drc)
quantile(Y.01.drc, c(.025, .975))

# Y(M(0), T=0)
Y.00.drc <- apply(predict.sp.0$draws.predict[[1]], 2, mean)
mean(Y.00.drc)
quantile(Y.00.drc, c(.025, .975))

# ACME
mean(Y.11.drc-Y.01.drc)
quantile((Y.11.drc-Y.01.drc), c(.025, .975))

mean(Y.10.drc-Y.00.drc)
quantile((Y.10.drc-Y.00.drc), c(.025, .975))

mean(((Y.11.drc-Y.01.drc)+(Y.10.drc-Y.00.drc))/2)
quantile((((Y.11.drc-Y.01.drc)+(Y.10.drc-Y.00.drc))/2), c(.025, .975))

# ANDE
mean(Y.11.drc-Y.10.drc)
quantile((Y.11.drc-Y.10.drc), c(.025, .975))

mean(Y.01.drc-Y.00.drc)
quantile((Y.01.drc-Y.00.drc), c(.025, .975))

mean(((Y.11.drc-Y.10.drc)+(Y.01.drc-Y.00.drc))/2)
quantile((((Y.11.drc-Y.10.drc)+(Y.01.drc-Y.00.drc))/2), c(.025, .975))

#########################################################################################
### Sri Lanka

alpha(cbind(sri$E2a, sri$E2b, sri$E2c, sri$E2d, sri$E2e, sri$E2f, sri$E2g, sri$E2h, sri$E2i, sri$E2j ))

# Mediator Equation
sri.sv.ptg <- ictreg.joint(Y ~ female + age.z + edu.z + income.z + hh_size.z +  killed + trauma + eastern, 
                         J=3, 
                         data=sri[is.na(sri$prior)==F, ], 
                         treat="treatment",
                         outcome="PTG",
                         outcome.reg="linear",
                         constrained=TRUE,
                         maxIter = 20000)
summary(sri.sv.ptg)

# Outcome Equation
sri.ptg.sp <- ictreg.joint(Y ~ female + age.z + edu.z + income.z + hh_size.z +  killed + trauma + eastern + prior + PTG, 
                         J=3, 
                         data=sri, 
                         treat="treatment",
                         outcome="soc_part.2",
                         outcome.reg="linear",
                         constrained=TRUE,
                         maxIter = 30000)
summary(sri.ptg.sp)

# Predicted Values
predict.ptg <- predict.ictreg.joint(sri.sv.ptg, sensitive.value="both", return.draws=T)

X.1 <- as.data.frame(sri.ptg.sp$x)
X.1$PTG <- predict.ptg[[1]]$Z1[,1]

X.0 <- as.data.frame(sri.ptg.sp$x)
X.0$PTG <- predict.ptg[[1]]$Z0[,1]

predict.sp.1 <- predict.ictreg.joint(sri.ptg.sp, newdata=X.1, sensitive.value="both", return.draws=T)
predict.sp.0 <- predict.ictreg.joint(sri.ptg.sp, newdata=X.0, sensitive.value="both", return.draws=T)

# Y(M=1, T=1)
Y.11.sri <- mean(predict.sp.1[[1]]$Z1[,1])

# Y(M=1, T=0)
Y.10.sri <- mean(predict.sp.1[[1]]$Z0[,1])

# Y(M=0, T=1)
Y.01.sri <- mean(predict.sp.0[[1]]$Z1[,1])

# Y(M=0, T=0)
Y.00.sri <- mean(predict.sp.0[[1]]$Z0[,1])

# ACME
Y.11.sri-Y.01.sri
Y.10.sri-Y.00.sri
((Y.11.sri-Y.01.sri)+(Y.10.sri-Y.00.sri))/2

# ANDE
Y.11.sri-Y.10.sri
Y.01.sri-Y.00.sri
((Y.11.sri-Y.10.sri)+(Y.01.sri-Y.00.sri))/2

# With Uncertainties
# Y(M(1), T=1)
Y.11.sri <- apply(predict.sp.1$draws.predict[[2]], 2, mean)
mean(Y.11.sri)
quantile(Y.11.sri, c(.025, .975))

# Y(M(1), T=0)
Y.10.sri <-  apply(predict.sp.1$draws.predict[[1]], 2, mean)
mean(Y.10.sri)
quantile(Y.10.sri, c(.025, .975))

# Y(M(0), T=1)
Y.01.sri <- apply(predict.sp.0$draws.predict[[2]], 2, mean)
mean(Y.01.sri)
quantile(Y.01.sri, c(.025, .975))

# Y(M(0), T=0)
Y.00.sri <- apply(predict.sp.0$draws.predict[[1]], 2, mean)
mean(Y.00.sri)
quantile(Y.00.sri, c(.025, .975))

# ACME
mean(Y.11.sri-Y.01.sri)
quantile((Y.11.sri-Y.01.sri), c(.025, .975))

mean(Y.10.sri-Y.00.sri)
quantile((Y.10.sri-Y.00.sri), c(.025, .975))

mean(((Y.11.sri-Y.01.sri)+(Y.10.sri-Y.00.sri))/2)
quantile((((Y.11.sri-Y.01.sri)+(Y.10.sri-Y.00.sri))/2), c(.025, .975))

# ANDE
mean(Y.11.sri-Y.10.sri)
quantile((Y.11.sri-Y.10.sri), c(.025, .975))

mean(Y.01.sri-Y.00.sri)
quantile((Y.01.sri-Y.00.sri), c(.025, .975))

mean(((Y.11.sri-Y.10.sri)+(Y.01.sri-Y.00.sri))/2)
quantile((((Y.11.sri-Y.10.sri)+(Y.01.sri-Y.00.sri))/2), c(.025, .975))

#########################################################################################
